# load dependencies
library(rgdal)
Loading required package: sp
rgdal: version: 1.4-8, (SVN revision 845)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE
Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
Linking to sp version: 1.3-2
library(geojsonio)
package ‘geojsonio’ was built under R version 3.6.2
Attaching package: ‘geojsonio’
The following object is masked from ‘package:base’:
pretty
library(ggplot2)
Registered S3 method overwritten by 'dplyr':
method from
print.rowwise_df
library(leaflet)
library(plotly)
package ‘plotly’ was built under R version 3.6.2Registered S3 method overwritten by 'data.table':
method from
print.data.table
Attaching package: ‘plotly’
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
library(RColorBrewer)
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(broom)
library(mapproj)
Loading required package: maps
library(reshape2)
package ‘reshape2’ was built under R version 3.6.2
library(digest)
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:base’:
date
library(spdplyr)
library(rmapshaper)
package ‘rmapshaper’ was built under R version 3.6.2Registered S3 method overwritten by 'geojsonlint':
method from
print.location dplyr
library(raster)
package ‘raster’ was built under R version 3.6.2
Attaching package: ‘raster’
The following object is masked from ‘package:dplyr’:
select
The following object is masked from ‘package:plotly’:
select
library(mapview)
package ‘mapview’ was built under R version 3.6.2
library(sf)
package ‘sf’ was built under R version 3.6.2Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
# try this with: https://www.r-graph-gallery.com/168-load-a-shape-file-into-r.html
zipcode_boundaries <- readOGR("../ZIP_CODE_040114/ZIP_CODE_040114.shp",
layer = "ZIP_CODE_040114", verbose = FALSE)
# try different table of zipcodes:
#alt_zipcode_boundaries <- readOGR("../alt_nyczip/nyc_zip_code_tabulation_areas_polygons.shp",
# layer = "nyc_zip_code_tabulation_areas_polygons", verbose = FALSE)
# fortify the data to obtain dataframe for use by ggplot2
zipcode_boundaries_fortified <- tidy(zipcode_boundaries, region = "ZIPCODE")
# plot it
ggplot() +
geom_polygon(data = zipcode_boundaries_fortified,
aes(x = long, y = lat, group = group),
fill = "#69b3a2", color="white") +
theme_void()
# try using info from: https://nceas.github.io/oss-lessons/spatial-data-gis-law/3-mon-intro-gis-in-r.html
class(zipcode_boundaries)
extent(zipcode_boundaries)
plot(zipcode_boundaries, main = "NYC zipcode boundaries")
# convert spatial data frame to geojson
# nyc_json <- geojson_json(zipcode_boundaries)
# save to local file
# geojson_write(nyc_json, file = "./nyc.geojson")
alt_nyc_json <- geojson_json(alt_zipcode_boundaries)
geojson_write(alt_nyc_json, file = "./alt_nyc.geojson")
# create leaflet object
pal <- colorNumeric("viridis", NULL)
leaflet(alt_zipcode_boundaries) %>%
# addProviderTiles(providers$OpenStreetMap) %>%
addTiles() %>% #add basemap
# addPolygons(stroke = FALSE, smoothFactor = 0.3, fillOpacity = 1,
# fillColor = ~pal(borough),
# label = borough) %>%
addPolylines(color = "#444444",
weight = 1,
smoothFactor = 0.5,
opacity = 1.0)
### run it all again but this time with SF (faster?)
zipcode_sf <- st_read("alt_nyczip/nyc_zip_code_tabulation_areas_polygons.shp")
Reading layer `nyc_zip_code_tabulation_areas_polygons' from data source `/Users/changhoyoon/Documents/Chang_Ho_Work/Master/Coursework/BMI706 data visualisation/viz_proj/COVID-19_NYC_706/alt_nyczip/nyc_zip_code_tabulation_areas_polygons.shp' using driver `ESRI Shapefile'
Simple feature collection with 262 features and 12 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: -74.25576 ymin: 40.49584 xmax: -73.6996 ymax: 40.91517
CRS: 4326
altzip_sf <- st_read("../ZIP_CODE_040114/ZIP_CODE_040114.shp")
Reading layer `ZIP_CODE_040114' from data source `/Users/changhoyoon/Documents/Chang_Ho_Work/Master/Coursework/BMI706 data visualisation/viz_proj/ZIP_CODE_040114/ZIP_CODE_040114.shp' using driver `ESRI Shapefile'
Simple feature collection with 263 features and 12 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 913129 ymin: 120020.9 xmax: 1067494 ymax: 272710.9
CRS: 2263
ggplot() +
geom_sf(data=zipcode_sf, aes(fill = borough))

# mapview function
# mapview(zipcode_sf)
mapview(altzip_sf)
# try leaflet with this
pal <- colorFactor(rainbow(5), zipcode_sf$borough)
leaflet(zipcode_sf) %>%
addProviderTiles(providers$CartoDB.Positron) %>% #add simpler basemap
addPolygons(stroke = FALSE,
smoothFactor = 0.3,
fillOpacity = 0.2,
highlight = highlightOptions(color = "white", weight = 10, bringToFront = TRUE)) %>%
addPolylines(color = "#444444",
weight = 1,
smoothFactor = 0.5,
opacity = 1.0,
highlightOptions = highlightOptions(color = "white",
weight = 5,
bringToFront = F,
opacity = 1))
### combine zipcode geometric data with cleaned dataframe
#cleaned_sf <- geojsonsf::geojson_sf("combined_covid_data.geojson")
#colnames(zipcode_sf)[which(names(zipcode_sf) == "postalcode")] <- "ZIPCODE"
# make certain of zipcode overlap between the two dataframes
setdiff(cleaned_sf$ZIPCODE, altzip_sf$ZIPCODE)
character(0)
# altzip only need ZIPCODE and geometry columns:
altzip_sf_compact <- subset(altzip_sf, select = c('ZIPCODE', 'geometry'))
combined_sf <- merge(cleaned_sf, altzip_sf_compact, by = 'ZIPCODE')
Error in merge.sf(cleaned_sf, altzip_sf_compact, by = "ZIPCODE") :
merge on two sf objects not supported
#
cleaned_sf_compact <- within(cleaned_sf, rm('geometry'))
cleaned_sf_compact <- cleaned_sf_compact %>% st_set_geometry(NULL)
class(cleaned_sf_compact)
[1] "data.frame"
#combined_sf <- merge(cleaned_sf, altzip_sf_compact, by = 'ZIPCODE')
#combined_sf <- merge(altzip_sf_compact, cleaned_sf_compact, by = 'ZIPCODE')
class(combined_sf)
[1] "sf" "data.frame"
mapview(combined_sf)
# let's try leaflet with this new combined_sf:
pal <- colorFactor(rainbow(5), combined_sf$PO_NAME)
leaflet(combined_sf) %>%
addProviderTiles(providers$CartoDB.Positron) %>% #add simpler basemap
addPolygons(stroke = FALSE,
smoothFactor = 0.3,
fillOpacity = 0.2,
highlight = highlightOptions(color = "white", weight = 10, bringToFront = TRUE))
sf layer is not long-lat datasf layer has inconsistent datum (+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs ).
Need '+proj=longlat +datum=WGS84'
# addPolylines(color = "#444444",
# weight = 1,
# smoothFactor = 0.5,
# opacity = 1.0,
# highlightOptions = highlightOptions(color = "white",
# weight = 5,
# bringToFront = F,
# opacity = 1))
LS0tCnRpdGxlOiBjb3ZpZDE5IFIgU2hpbnkgbm90ZWJvb2sKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CiMgbG9hZCBkZXBlbmRlbmNpZXMKbGlicmFyeShyZ2RhbCkKbGlicmFyeShnZW9qc29uaW8pCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShsZWFmbGV0KQpsaWJyYXJ5KHBsb3RseSkKbGlicmFyeShSQ29sb3JCcmV3ZXIpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkobWFwcHJvaikKbGlicmFyeShyZXNoYXBlMikKbGlicmFyeShkaWdlc3QpCmxpYnJhcnkobHVicmlkYXRlKQpsaWJyYXJ5KHNwZHBseXIpCmxpYnJhcnkocm1hcHNoYXBlcikKbGlicmFyeShyYXN0ZXIpCmxpYnJhcnkobWFwdmlldykKbGlicmFyeShzZikKYGBgCgoKYGBge3J9CiMgdHJ5IHRoaXMgd2l0aDogaHR0cHM6Ly93d3cuci1ncmFwaC1nYWxsZXJ5LmNvbS8xNjgtbG9hZC1hLXNoYXBlLWZpbGUtaW50by1yLmh0bWwKCnppcGNvZGVfYm91bmRhcmllcyA8LSByZWFkT0dSKCIuLi9aSVBfQ09ERV8wNDAxMTQvWklQX0NPREVfMDQwMTE0LnNocCIsCiAgbGF5ZXIgPSAiWklQX0NPREVfMDQwMTE0IiwgdmVyYm9zZSA9IEZBTFNFKQoKIyB0cnkgZGlmZmVyZW50IHRhYmxlIG9mIHppcGNvZGVzOgojYWx0X3ppcGNvZGVfYm91bmRhcmllcyA8LSByZWFkT0dSKCIuLi9hbHRfbnljemlwL255Y196aXBfY29kZV90YWJ1bGF0aW9uX2FyZWFzX3BvbHlnb25zLnNocCIsCiMgIGxheWVyID0gIm55Y196aXBfY29kZV90YWJ1bGF0aW9uX2FyZWFzX3BvbHlnb25zIiwgdmVyYm9zZSA9IEZBTFNFKQoKIyBmb3J0aWZ5IHRoZSBkYXRhIHRvIG9idGFpbiBkYXRhZnJhbWUgZm9yIHVzZSBieSBnZ3Bsb3QyCnppcGNvZGVfYm91bmRhcmllc19mb3J0aWZpZWQgPC0gdGlkeSh6aXBjb2RlX2JvdW5kYXJpZXMsIHJlZ2lvbiA9ICJaSVBDT0RFIikKCiMgcGxvdCBpdApnZ3Bsb3QoKSArCiAgZ2VvbV9wb2x5Z29uKGRhdGEgPSB6aXBjb2RlX2JvdW5kYXJpZXNfZm9ydGlmaWVkLAogICAgICAgICAgICAgICBhZXMoeCA9IGxvbmcsIHkgPSBsYXQsIGdyb3VwID0gZ3JvdXApLAogICAgICAgICAgICAgICBmaWxsID0gIiM2OWIzYTIiLCBjb2xvcj0id2hpdGUiKSArCiAgdGhlbWVfdm9pZCgpCmBgYAoKYGBge3J9CiMgdHJ5IHVzaW5nIGluZm8gZnJvbTogaHR0cHM6Ly9uY2Vhcy5naXRodWIuaW8vb3NzLWxlc3NvbnMvc3BhdGlhbC1kYXRhLWdpcy1sYXcvMy1tb24taW50cm8tZ2lzLWluLXIuaHRtbAoKY2xhc3MoemlwY29kZV9ib3VuZGFyaWVzKQoKZXh0ZW50KHppcGNvZGVfYm91bmRhcmllcykKCnBsb3QoemlwY29kZV9ib3VuZGFyaWVzLCBtYWluID0gIk5ZQyB6aXBjb2RlIGJvdW5kYXJpZXMiKQpgYGAKCgpgYGB7cn0KIyBjb252ZXJ0IHNwYXRpYWwgZGF0YSBmcmFtZSB0byBnZW9qc29uCiMgbnljX2pzb24gPC0gZ2VvanNvbl9qc29uKHppcGNvZGVfYm91bmRhcmllcykKIyBzYXZlIHRvIGxvY2FsIGZpbGUKIyBnZW9qc29uX3dyaXRlKG55Y19qc29uLCBmaWxlID0gIi4vbnljLmdlb2pzb24iKQoKYWx0X255Y19qc29uIDwtIGdlb2pzb25fanNvbihhbHRfemlwY29kZV9ib3VuZGFyaWVzKQpnZW9qc29uX3dyaXRlKGFsdF9ueWNfanNvbiwgZmlsZSA9ICIuL2FsdF9ueWMuZ2VvanNvbiIpCmBgYAoKYGBge3J9CiMgY3JlYXRlIGxlYWZsZXQgb2JqZWN0CnBhbCA8LSBjb2xvck51bWVyaWMoInZpcmlkaXMiLCBOVUxMKQoKbGVhZmxldChhbHRfemlwY29kZV9ib3VuZGFyaWVzKSAlPiUKIyAgYWRkUHJvdmlkZXJUaWxlcyhwcm92aWRlcnMkT3BlblN0cmVldE1hcCkgJT4lIAogIGFkZFRpbGVzKCkgJT4lICNhZGQgYmFzZW1hcAojICBhZGRQb2x5Z29ucyhzdHJva2UgPSBGQUxTRSwgc21vb3RoRmFjdG9yID0gMC4zLCBmaWxsT3BhY2l0eSA9IDEsCiMgICAgZmlsbENvbG9yID0gfnBhbChib3JvdWdoKSwKIyAgICBsYWJlbCA9IGJvcm91Z2gpICU+JQogIGFkZFBvbHlsaW5lcyhjb2xvciA9ICIjNDQ0NDQ0IiwgCiAgICAgICAgICAgICAgIHdlaWdodCA9IDEsIAogICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjUsIAogICAgICAgICAgICAgICBvcGFjaXR5ID0gMS4wKQpgYGAKCgpgYGB7cn0KIyMjIHJ1biBpdCBhbGwgYWdhaW4gYnV0IHRoaXMgdGltZSB3aXRoIFNGIChmYXN0ZXI/KQoKemlwY29kZV9zZiA8LSBzdF9yZWFkKCJhbHRfbnljemlwL255Y196aXBfY29kZV90YWJ1bGF0aW9uX2FyZWFzX3BvbHlnb25zLnNocCIpCmFsdHppcF9zZiA8LSBzdF9yZWFkKCIuLi9aSVBfQ09ERV8wNDAxMTQvWklQX0NPREVfMDQwMTE0LnNocCIpCgpnZ3Bsb3QoKSArCiAgZ2VvbV9zZihkYXRhPXppcGNvZGVfc2YsIGFlcyhmaWxsID0gYm9yb3VnaCkpIApgYGAKCmBgYHtyfQojIG1hcHZpZXcgZnVuY3Rpb24KIyBtYXB2aWV3KHppcGNvZGVfc2YpCm1hcHZpZXcoYWx0emlwX3NmKQpgYGAKCmBgYHtyfQojIHRyeSBsZWFmbGV0IHdpdGggdGhpcwpwYWwgPC0gY29sb3JGYWN0b3IocmFpbmJvdyg1KSwgemlwY29kZV9zZiRib3JvdWdoKQoKbGVhZmxldCh6aXBjb2RlX3NmKSAlPiUKICBhZGRQcm92aWRlclRpbGVzKHByb3ZpZGVycyRDYXJ0b0RCLlBvc2l0cm9uKSAlPiUgICNhZGQgc2ltcGxlciBiYXNlbWFwIAogIGFkZFBvbHlnb25zKHN0cm9rZSA9IEZBTFNFLCAKICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjMsIAogICAgICAgICAgICAgIGZpbGxPcGFjaXR5ID0gMC4yLAogICAgaGlnaGxpZ2h0ID0gaGlnaGxpZ2h0T3B0aW9ucyhjb2xvciA9ICJ3aGl0ZSIsIHdlaWdodCA9IDEwLCBicmluZ1RvRnJvbnQgPSBUUlVFKSkgICU+JQogIGFkZFBvbHlsaW5lcyhjb2xvciA9ICIjNDQ0NDQ0IiwgCiAgICAgICAgICAgICAgIHdlaWdodCA9IDEsIAogICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjUsIAogICAgICAgICAgICAgICBvcGFjaXR5ID0gMS4wLCAKICAgICAgICAgICAgICAgaGlnaGxpZ2h0T3B0aW9ucyA9IGhpZ2hsaWdodE9wdGlvbnMoY29sb3IgPSAid2hpdGUiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0ID0gNSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJyaW5nVG9Gcm9udCA9IEYsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBvcGFjaXR5ID0gMSkpCmBgYAoKYGBge3J9CiMjIyBjb21iaW5lIHppcGNvZGUgZ2VvbWV0cmljIGRhdGEgd2l0aCBjbGVhbmVkIGRhdGFmcmFtZQoKI2NsZWFuZWRfc2YgPC0gZ2VvanNvbnNmOjpnZW9qc29uX3NmKCJjb21iaW5lZF9jb3ZpZF9kYXRhLmdlb2pzb24iKQoKI2NvbG5hbWVzKHppcGNvZGVfc2YpW3doaWNoKG5hbWVzKHppcGNvZGVfc2YpID09ICJwb3N0YWxjb2RlIildIDwtICJaSVBDT0RFIgoKIyBtYWtlIGNlcnRhaW4gb2YgemlwY29kZSBvdmVybGFwIGJldHdlZW4gdGhlIHR3byBkYXRhZnJhbWVzCnNldGRpZmYoY2xlYW5lZF9zZiRaSVBDT0RFLCBhbHR6aXBfc2YkWklQQ09ERSkKCiMgYWx0emlwIG9ubHkgbmVlZCBaSVBDT0RFIGFuZCBnZW9tZXRyeSBjb2x1bW5zOgphbHR6aXBfc2ZfY29tcGFjdCA8LSBzdWJzZXQoYWx0emlwX3NmLCBzZWxlY3QgPSBjKCdaSVBDT0RFJywgJ2dlb21ldHJ5JykpCgpgYGAKCgpgYGB7cn0KIyAKY2xlYW5lZF9zZl9jb21wYWN0IDwtIHdpdGhpbihjbGVhbmVkX3NmLCBybSgnZ2VvbWV0cnknKSkKY2xlYW5lZF9zZl9jb21wYWN0IDwtIGNsZWFuZWRfc2ZfY29tcGFjdCAlPiUgc3Rfc2V0X2dlb21ldHJ5KE5VTEwpCmNsYXNzKGNsZWFuZWRfc2ZfY29tcGFjdCkKI2NvbWJpbmVkX3NmIDwtIG1lcmdlKGNsZWFuZWRfc2YsIGFsdHppcF9zZl9jb21wYWN0LCBieSA9ICdaSVBDT0RFJykKYGBgCgpgYGB7cn0KI2NvbWJpbmVkX3NmIDwtIG1lcmdlKGFsdHppcF9zZl9jb21wYWN0LCBjbGVhbmVkX3NmX2NvbXBhY3QsIGJ5ID0gJ1pJUENPREUnKQoKY2xhc3MoY29tYmluZWRfc2YpCgptYXB2aWV3KGNvbWJpbmVkX3NmKQpgYGAKCmBgYHtyfQojIGxldCdzIHRyeSBsZWFmbGV0IHdpdGggdGhpcyBuZXcgY29tYmluZWRfc2Y6CgpwYWwgPC0gY29sb3JGYWN0b3IocmFpbmJvdyg1KSwgY29tYmluZWRfc2YkUE9fTkFNRSkKCmxlYWZsZXQoY29tYmluZWRfc2YpICU+JQogIGFkZFByb3ZpZGVyVGlsZXMocHJvdmlkZXJzJENhcnRvREIuUG9zaXRyb24pICU+JSAgI2FkZCBzaW1wbGVyIGJhc2VtYXAgCiAgYWRkUG9seWdvbnMoc3Ryb2tlID0gRkFMU0UsIAogICAgICAgICAgICAgIHNtb290aEZhY3RvciA9IDAuMywgCiAgICAgICAgICAgICAgZmlsbE9wYWNpdHkgPSAwLjIsCiAgICBoaWdobGlnaHQgPSBoaWdobGlnaHRPcHRpb25zKGNvbG9yID0gIndoaXRlIiwgd2VpZ2h0ID0gMTAsIGJyaW5nVG9Gcm9udCA9IFRSVUUpKQojICBhZGRQb2x5bGluZXMoY29sb3IgPSAiIzQ0NDQ0NCIsIAojICAgICAgICAgICAgICAgd2VpZ2h0ID0gMSwgCiMgICAgICAgICAgICAgICBzbW9vdGhGYWN0b3IgPSAwLjUsIAojICAgICAgICAgICAgICAgb3BhY2l0eSA9IDEuMCwgCiMgICAgICAgICAgICAgICBoaWdobGlnaHRPcHRpb25zID0gaGlnaGxpZ2h0T3B0aW9ucyhjb2xvciA9ICJ3aGl0ZSIsIAojICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0ID0gNSwgCiMgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBicmluZ1RvRnJvbnQgPSBGLCAKIyAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG9wYWNpdHkgPSAxKSkKYGBgCgo=